SUBROUTINE rmultinomial(u, n, prob, out)
  USE Commonvars
  IMPLICIT NONE
  REAL(8), INTENT(IN)  :: u
  INTEGER, INTENT(IN)  :: n
  REAL(8), INTENT(IN)  :: prob(n)
  INTEGER, INTENT(OUT) :: out
  REAL(8)              :: x
  INTEGER              :: i

  IF(u > 1d0 .OR. u < 0d0) THEN
     WRITE(*,*) 'u must be between 0 and 1 (rmultinomial)'
     write(*,*) u
     IF(istop == 1) STOP 
  END IF

  IF(SUM(prob) < 1d0 - 1d-2 .OR. SUM(prob) > 1d0 + 1d-2) THEN
     WRITE(*,*) 'prob must sum 1. sum(prob) is ', sum(prob),&
          & SUM(prob) - 1d0
     WRITE(*,*) prob
     IF(istop == 1) STOP
  END IF

  x = 0d0
  DO i = 1, n
     x = x + prob(i)
     IF(u < x) THEN
        out = i
        RETURN
     END IF
  END DO

  out = n

END SUBROUTINE rmultinomial

!!$PROGRAM main
!!$  INTEGER :: x
!!$
!!$  CALL rmultinomial(.51d0, 4, (/.25d0, .25d0, .25d0, .25d0/), x)
!!$  WRITE(*,*) x
!!$END PROGRAM main
